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Abstract 

The volume penalty method provides a simple, efficient approach for solving 
the incompressible Navier-Stokes equations in domains with boundaries or 
in the presence of moving objects. Despite the simplicity, the method is 
typically limited to first order spatial accuracy. We demonstrate that one 
may achieve high order accuracy by introducing an active penalty term. One 
key difference from other works is that we use a sharp, unregularized mask 
function. We discuss how to construct the active penalty term, and provide 
numerical examples, in dimensions one and two. We demonstrate second and 
third order convergence for the heat equation, and second order convergence 
for the Navier-Stokes equations. In addition, we show that modifying the 
penalty term does not significantly alter the time step restriction from that 
of the conventional penalty method. 

Keywords: Active penalty method, Sharp mask function, Immersed 
boundary, Incompressible flow, Navier-Stokes, Heat equation 



1. Introduction 

There are many popular methods for numerically solving the incompress- 
ible Navier-Stokes equations in complex geometries. For instance, the im- 
mersed boundary method |20j, the immersed interface method [18] and the 
ghost fluid method pj] are popular since they allow one to use a regular grid 
with an immersed domain boundary. In such cases, the regular grid alleviates 
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many of the numerical difficulties introduced by curved or moving bound- 
aries. In this paper, we focus on the volume penalty method [TJ EJ El [7J HE] , 
which loosely fits into the same class of methods. 

As a result of their simplicity, penalty methods have been used in a wide 
variety of problems including electromagnetism, magnetohydrodynamics [Tj5] , 
shape optimization [9], fluid-solid interaction problems [TO], [T3] and even sim- 
pler problems such as the heat equation or Poisson equation [2T]. In the 
context of fluids, they provide a simple means for solving the incompressible 
Navier-Stokes equations in domains with boundaries. The approach relies 
on replacing the often difficult to implement Dirichlet fluid boundary condi- 
tions, with a simpler to implement volumetric forcing term in the advection 
equation. 

Despite the simplicity, the penalty method suffers from i) poor conver- 
gence in the penalty parameter, thereby restricting the accuracy of numerical 
methods and, ii) a lack of regularity in the velocity field which reduces the 
advantages of spectral methods. For example, solutions to the penalized 
equations have a discontinuous second derivative which limits the decay rate 
of the Fourier coefficients, as well as the ability to spectrally compute deriva- 
tives. Despite the lack of smoothness, stable and low order spectral methods 
have been successfully used to solve the penalized fluid equations [151 E] • 

The focus of our paper is to introduce a systematic method for improving 
the accuracy of penalty methods. Current methods which improve accuracy 
rely on introducing a subgrid numerical construct in the vicinity of the do- 
main boundary [22| [23] . Such approaches, however, are restrictive if one 
wishes to eventually use high order Fourier methods. One distinct differ- 
ence with our approach is that we alter the equations at the continuous level 
to improve the analytic convergence rate of the penalized problem to the 
original unpenalized problem. The improved analytic convergence rate then 
allows for higher order numerical schemes. 

We first introduce the original volume penalty method, followed by an 
introduction to the improved active penalty method. We then explicitly show 
how to analytically construct the new penalty term. Following the con- 
struction, we then examine a model equation to demonstrate how the active 
penalization improves the convergence rate for the Poisson equation. 

After discussing the improved convergence, we focus on numerical details. 
First, we examine the stability of the new active penalty term, and show that 
it does not introduce additional numerical stiffness. We then provide numer- 
ical examples for the heat equation, in dimensions one and two, showing 
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second and third order schemes. Lastly, we outline how to handle the di- 
vergence constraint for the Navier- Stokes equations and provide numerical 
examples showing second order convergence (in L°°) in the velocity field and 
first order in the pressure. 

2. Navier-Stokes and volume penalty equations 

The aim of our work is to examine the behavior of a fluid in the vicinity of 
a solid or a porous medium. For instance, two examples include the motion 
of a fluid in a bounded domain with hard walls, or the motion around an 
immersed solid body such as the one shown in figure [TJ In our case, we 
consider dimensions D = 2, 3 and let fl p C K, D denote the physical fluid 
domain. For our purposes, Q p is an open set with C 2 boundary Y = dVL p . 

2.1. Incompressible Navier-Stokes equations 

Through the conservation of mass and momentum, the incompressible 
Navier-Stokes equations govern the flow of an incompressible fluid for x G Q p 

d t u + u • Vu = -Vp + fiAu + f (1) 
V-u = 0. (2) 

Here u(x, t) is the velocity vector field, p(x, t) is the pressure, \x > is the 
kinetic viscosity, and f(x, t) is an external forcing such as gravity. 

To supplement the bulk equations (JT])-(j2]), the fluid velocity also satisfies 
prescribed boundary conditions 

u = g for x G T (3) 
J g • n dA = (4) 

Here n is an outward pointing normal, while equation Q represents a con- 
sistency condition on the boundary data. Although we allow g = g(x, t) 
to be a function of both space and time, the case of g = represents the 
practical condition of a no-slip and no-flux boundary condition. Together, 
equations Q-pl) with boundary data ^ describe the evolution of an initial, 
divergence-free velocity field u(x, 0) = u (x). 
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2.2. Volume penalty equations 

Domains with curved boundaries T present several challenges to the nu- 
merical solution of equations Q-([2]). For example, curved boundaries or 
immersed objects limit the use of Fourier methods since solutions are not 
periodic, or easily extended to periodic functions. One simple solution to 
handle complicated or moving boundaries is through the use of a volume 
penalty term in the Navier-Stokes equations. In such a case, one removes 
the Navier-Stokes boundary condition, and instead adds a drag term to the 
momentum equation. 

To introduce the penalized equation, we first denote Q s as the solid do- 
main of the obstacle or wall. Here the obstacle region Q s = Q s is a closed 
set which shares the same boundary as the fluid, dQ s = T. The penalized 
equations are then defined on a computational domain Q which is the union 
of the physical and solid domains Q = Q p [J Q s . In our case we take Q to be 
a rectangular domain with periodic boundary conditions, i.e. Q = T D where 
T D is the D-dimensional torus. 




Figure 1: Illustration of the physical fluid (O p ) and solid obstacle (fl s ) do- 
mains. 

For a stationary obstacle with a g = boundary condition, the volume 
penalty equations (see [31IHE]) are 

d t u v + u, n ■ Vu r/ = — Vprj + fiAu ri + f — ri^Xs u r/ x G Vt (5) 
V-u, = 0. (6) 

Here rj is a small parameter, and Xs( x ) is the characteristic function on Q s , 
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namely 



Xs(x)= J° forxefi\fi s (7) 
x 1 1 for x£fl s 

In the limit 77 — > 0, the drag term in equation ^ becomes large and tends 
to slow the fluid inside Q s . Rigorous convergence results by Angot et al. [3], 
and Carbou and Fabrie [S] show that the penalized velocity u,, converges to 
the solution of the Navier-Stokes equations u with an error rate of Ofa 1 / 2 ) 
in the L 2 (Q p ) norm. 

2.3. Improved volume penalty equations 

Although the volume penalty equations do converge to Navier-Stokes as 
7] — > 0, the convergence rate is slow and therefore may limit the accuracy 
of resulting numerical schemes. For example, let u^^™ denote a numerical 
solution for the penalized equations. One is then interested in quantifying 
the numerical error for u^ nnm compared to u, the solution to the original 
Navier-Stokes problem ([TJ-Q. Using the triangle inequality^ the error can 
be controlled by 

I | u ~~ u ri,num\ h < | | u ~~ u r;| I2 + I | u ?7 ~~ u rj,num\ \%- (8) 

Rigorous convergence results then bound the first term as ||u — u^l^ ~ ?? 1 ^ 2 , 
while ||u^ — Ur/ : ?T.«m, 1 12 depends on the numerical details and order of the 
scheme. Finally, we note that the addition of the penalty term introduces 
time scales of 0(rf) and length scales of O^ 1 / 2 ) into the solution u,,. To 
appropriately resolve the boundary layers in the penalty equations (J5])-(|6|, 
one then has a grid spacing of Ax ~ ^ff] leading to a first order bound 

||u-u nirium || 2 < O(Ax). (9) 

In light of the above observations, a high order penalty method must either 
increase the boundary layer width 0(y/rj), or improve the analytic conver- 
gence rate in the penalty parameter. We adopt the second approach, and 
outline how equation (JsJ) can be modified to better approximate the original 
Navier-Stokes problem ([T])-([2]). Furthermore, we note that when modify- 
ing the penalty term, it is important to avoid the introduction of additional 



L Here || • H2 is any appropriate numerical L 2 (f2 p ) norm. 
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length or time scales which would hinder the development of high order nu- 
merical schemes. To improve the penalized equations, we exploit the fact that 
u satisfies the boundary conditions on T, and does not represent a physical 
flow inside fl s . Specifically, we modify the penalty term so that the flow 
tracks an extension function g defined on x G Q s . In such a case, the volume 
penalty equations take the form 

d t u v + u,, ■ Vu, = -Vp v + fiAu n + f - rj^Xs (u v - g) xed (10) 
V-u„ = xefip. (11) 

At this point, we only specify the divergence constraint within the physi- 
cal domain and defer a more detailed description of the divergence constraint 
inside Q s for section [7j The idea is to choose g to reduce the artificial fluid 
boundary layer generated by the penalized equations in the vicinity of T. 
Specifically, the function g should be a smooth, at least C 1 , extension of 
the prescribed boundary conditions. The extension is constructed for each 
component of g, and each component of g should be chosen to match k 
derivatives of u in the direction normal to T. As a result, we prescribe the 
following general properties for g 

PI. g is an extension of the prescribed boundary values: g = g for x G T. 
P2. g has the same normal slope as u: (n ■ V)?ij = (n ■ V)(?j. Here Ui and 

gi for i = 1 . . . D are the components of u and g 
P3. For higher derivatives, (n ■ V) k Ui\r p = (n • V) fc ^. 

Since derivatives of u may be discontinuous across T, the notation T p denotes 
the limit of the derivative from the physical domain Q p . 

3. Constructing the extension g 

In this section we discuss one possible construction for the extension 
function g. The construction procedure is identical for each component gi of 
g for i = 1 . . . D. 

In our approach, we assume the domain Q p has a smooth boundary, at 
least r G C 2 . As a result, we omit a class of physically important domains 
such as rectangles. The general idea is to match the normal derivatives of g 
to those of u on T. With the appropriate boundary derivatives, we then let 
g decay to some constant value G over a length scale I. In our construction, 
the maximum length scale I is bounded by the minimum radius of curvature 
of the interface. 



6 



Step 1. First introduce a family of smooth, one-dimensional basis functions 
Bj G C k with < j < k such that 

(i) The functions Bj form a basis for derivatives at x = 

d 1 , . I 1 for ?' = 2 , , 

— R = < J 12 

dx* n ' |0 for j ^ i y ' 

(ii) Each Bj(x) has support on [0,1]. Namely Bj(x) = for x < 
and x > 1. 

One can then use the functions Bj(x) to construct a C fc extension /(x) 
of a one-dimensional function /(x) on x > as 

/>) = /(O)Sq(x) + /'(0)5!(x) + . . . + /«(0)S fe (x). (13) 

Note that by construction, the function f(x) matches k derivatives at 
x = and vanishes for x > 1. The goal is now to modify the extension 
(13) to higher dimensions. 

Although there are many different choices for Bj(x), we now give an 
example of one such choice for matching k = 2 derivatives. We do this 
by constructing Bj(x) out of stretched copies of the smoothly decaying 
function 

h(x) = l el ~^ fOT 0<x<l^ 
for x > 1 



Using h(x), one can take the functions B , B l ,B 2 (figure [2]) as the 
weighted sums 

B (x) = 3 h(x) - 3 h(2x) + h(3x) (15) 
B 1 {x) = - h(x) - 4 h(2x) + ^ h(3x) (16) 

B 2 {x) = -- h{x) + h(2x) - - h(3x). (17) 

Step 2. Construct a coordinate system inside the obstacle. The coordinate 
system should be orthogonal at the boundary, and only needs to extend 
a distance / inside the domain Q s . 

To construct the coordinates, we follow a standard approach [12] shown 
in figure [3j Let £ G T denote the coordinates of the boundary T. Any 
point x within a distance / of the boundary may then be written as 

x = £ + S n(0- (18) 
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Figure 2: A plot of the ID basis functions B (x) (thick), B\(x) (dashed) and 
B 2 (x) (thin). 



Here s is the normal distance inside f2 s from the boundary. Within 
a small enough region, s < I, one may invert^] equation (18) to write 
£ = £(x) and s = s(x). 

Remark 1. In cases where a level set function 0(x) describes the 
boundary r = {x G H <K X ) = 0}? one m &y identify 

n = V0 (19) 
s = 0(x). (20) 

Here we have assumed |V0| = 1 and 0(x) > represents the domain 
Q s while 0(x) < corresponds to the domain Q p . 6 

Step 3. Construct the extension g using the functions Bj(x) and the coordi- 
nates (£, s). 

For brevity, we introduce notation for the normal derivatives at the 



2 The coordinates £(x) and s(x) are both at least C 1 functions 
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Figure 3: A local set of normal coordinates. Here $, G T is a point on the 
boundary, while s is the distance in the normal direction. A coordinate inside 
a neighborhood of fl s has the form x = £ + s n(£). 



boundary T. 

u n (£) = (n-V)u| r (21) 
u„„(£) = (n-V) 2 u| r (22) 

u nfc (£) = (n-V)Si| r (23) 

Again, we note that higher derivatives of u are discontinuous across 
the boundary. Therefore, u nk is evaluated as the limit approaching the 
boundary from the physical domain. The extension is then 

g(x) = (g(O-G) B (s r x )+l u n (£) B l (sT 1 ) 
+ I 2 u nn (|) B 2 (sl- 1 )+... 

+ l k u nk (£) B^sT^+G. (24) 

Note that g decays to G, i.e. g — > G, as s — > I. Therefore G can be 
any time-dependent constant vector, however, for numerical purposes 
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one should choose G close to the boundary average of g 



G = A' 1 j gdA (25) 
A = J dA (26) 
Remark 2. Since values of g inside fl s depend on derivatives of u on the 



boundary, the function g described in (24) depends linearly on u. 4k 



Remark 3. We can check that the construction (24) satisfies the properties 



(P1)-(P3). For x G T, we have s = 0, so that g = g, thereby satisfying 



(PI). To check higher derivatives, we first note that differentiating (18) with 
respect to s yields = n(£). As a result any function independent of s, 
y($,,s) = ?/(£), has the property that 

, % (a8) 

= 0. (29) 

Meanwhile, we also have 

(n-V^sZ- 1 )!^ = (—YB^sr 1 )^ (30) 



ds< 

= I'^ij (31) 

where 5ij is the Kronecker delta. Combining the two properties above, we 
have 

(n-V) l g| xer = u nj (0- (32) 
Therefore, we recover properties (P2)-(P3). 



4. A Model equation 



In this section we examine solutions to the steady-state heat equation 
to provide some explanation for how the extension function g improves the 
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analytic convergence rate of the penalized equations to the original problem. 
In particular, we seek to quantify the error that results from the additional 
penalty forcing. As a non-penalized problem, consider 

d xx v = x G [-1,0] (33) 

with boundary conditions: i>(— 1) = 1, v (0) = 0. The solution is then v(x) = 
—x for x G [—1,0]. 

We note that solving explicit examples of the steady-state equations do 
not give general sharp convergence estimates, however, they do provide a 
rigorous lower bound on the convergence rate of the penalized equation to 
the exact non-penalized equation. The equivalent one dimensional steady- 
state penalized problem is then 

d^u = r]~ 1 H(x)(u - g) xG[-l,oo) (34) 

with boundary conditions u(—l) = 1, u(oo) = 0. Here H{x) is the Heaviside 
function 

H(x) = f° ioix<0 (35) 
1 1 for x > 

We now examine the convergence of solutions u — > v in the limit r\ — > for 
different extensions g. 

Remark 4. As a result of the discontinuous Heaviside function H(x), the 



solution u to equation (34) has one continuous derivative (u G C 1 ). Higher 



derivatives are discontinuous across x — 0. 4k 

In light of remark we take g to have derivatives matching u from the 
physical domain x = 0~ and not x = + . 

Proposition 1. Suppose that g is a bounded C k+l function that matches k 
derivatives of u at x = 0". Namely 

1. ff(0)=0 

2. g'(0) = u'(0) 

3. (?( m )(0) = u^(Q-) = for 2 < m < k. 



11 



Then the solution u converges to v as 

max \u-v\ = 0(n {k+1)/2 ) (36) 
xe[-i,o} 

Proof 1. In the region — 1 < x < 0, u has the solution 

u = (1 + c) + cx (37) 

for some constant c. To construct the solution on x > 0, we note that one 
may write g as 

g(x) =cx + x k+1 R(x) (38) 

for some remainder function R(x), where in general R(0) ^ 0. By construc- 
tion (38) matches the first k derivatives of u at x = 0~ . In addition, we 
assume that g(x) and g'(x) are bounded, so that R(x) and R'(x) are also 
bounded functions. On x > 0, u then solves 

d xx u - r]- l u = -r]~ l {cx + x k+1 R(x)) (39) 

To obtain the correct scaling, we rescale x = rf^X to obtain 

d xxU - u = -c^X + V W 2 X k+1 R(riV 2 X) (40) 

The general solution is then 

{(l + ^ + c^X ifX<0 
H } \be- x + c V 1 / 2 X + r ] ( k+1 V 2 Q(X) if X > 

where we have excluded the term e x since it diverges as X — > oo. In addition, 
Q(X) is a particular solution (which stays bounded as X — )• oo ) to 

Qxx -Q = X k+1 R(^ 2 X) (42) 

For instance, one may write a particular solution as 

Q(X) = l -J X {e x -y-e- x+ y)y k+1 R{^ 2 y)Ay-Ae x (43) 
1 f°° 

A = -J e-yy k+l R( V ^ 2 y)dy (44) 
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Letting R m = max s \R(y)\, we also have the bound 

\Q(0)\ = \A\ <^J y k+1 e-ydy = Q . (45) 

The same type of argument holds for bounding |Q'(0)| < Q±. 

To solve for the unknown constants, c and b, we use the fact that u and 
u' are continuous across x = 0. We therefore obtain the two equations 

1 + c = b + 7] {k+1)/2 Q(0) (46) 

^/2 C = _ 6 + 77 V2 c + r ^+l)/2g/ (0) (47) 

Upon solving for b and c, the error between u and t> on the physical domain 
-1 < x < is 

max lit — u| = 1 + c (48) 

xel-i,o} 

= (Q(0) + Q'(0))?/ fc+1 >/ 2 (49) 

< (Qo + QiW k+1)/2 (50) 

= 0(77( fc+1 )/ 2 ) (51) 

Hence, for the model problem, matching k derivatives of g yields a conver- 
gence rate of (k + l)/2. In particular, when k = 0, we recover the known 
convergence rate 7/ 1//2 of the standard penalty method. 

Remark 5. Using higher derivatives in the construction of g which are taken 
as limits from the domain x \ + , does not yield the convergence rate stated 
in proposition ([!]). As an example, we take g + = u'(0)Bi(x) + u"(0 + )B2(x) 
where 

B 1 (x) = ^ e - x -4e- 2x + ^e- 3x (52) 
B 2 (x) = X -e~ x - e~ 2x + ^e~ 3x . (53) 



For such a g + , the solution u to problem (34) yields only a first order error 

11 . ^ 

max \u — v\ ~ — 77. (54) 
[-1,0]' 1 6 ' K ' 

In contrast, taking g~ = u'(0)Bi(x) + u"(0~)B2(x) yields a convergence rate 
in agreement with Q 

max \u — v\ ~ IIt? 3 ^ 2 . (55) 

[-1,0]' 
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5. Stability 



In this section we establish numerical stability criteria for the ID pe- 
nalized heat equation. To examine stability, we work with the domain 
Q p = (0, 7r), Q s = [vr,27r] and periodic boundary conditions. Moreover, we 
take g(0) = g(n) = to capture a u = boundary condition at the fluid- 
solid boundary. A simple Euler scheme matching one derivative of u at the 
interface is then 



" +1 = (/ + At A)u n - At-q^x, (u n -g n ) (56) 
g n = u£(7r)Bi(x-7r)-t£(27r)Bi(27r-aO (57) 



In general, adding derivatives of u to g(x) can reduce, by factors of h, the 
time step restriction for an explicit scheme. In the case at hand, however, 
the structure of g(x) results in the same time step restriction as the original 
volume penalty method, namely 

At <mm{0(h 2 ),0(ri)}. (58) 

Here h is either the grid spacing of a finite difference scheme, or alternatively 
h" 1 scales as the largest wavenumber in a Fourier method. 



We note that although (56) is a linear recursion relation, the right hand 



side is not a normal operator. As a result, a rigorous proof of (58) requires 



bounding the eigenvalues for the spatially discrete system (56). The analysis 



is further complicated by the fact that the operators (or matrices) on the 



right hand side of ( 56 ) do not commute 



In this section we establish the time step restriction (58). To do so, we 
first compute the eigenvalues for the penalty term using a finite difference 
scheme. We show that although the penalty term contains derivatives of u, 
the eigenvalues remain 0(r/ _1 ) and do not become 0(r/ _1 /i _1 ). 

Secondly, to show that the addition of the Laplacian does not alter the 



restriction (58), we numerically compute the eigenvalues for equation (56) 



using a finite difference scheme. 

5.1. Eigenvalues of the penalty term (Finite differences) 

In practice, one does not observe the time step restriction governed by the 



norm of g, but rather the larger bound in (58). Here we provide a stability 
criteria by analytically computing the penalty term eigenvalues for a finite 
difference scheme. 
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Let Xk = kh for < k < N — l with grid spacing h = 2tt/N. Furthermore, 
denote the discrete vector u = [u(xq) u{x\) . . . u(xn-i)] t ■ 

We are then interested in evaluating the eigenvalues of the penalty term 

Bu = Au (59) 
B = -7T 1 (I x -v 1 d? , -v 2 d!') (60) 

where B is the finite difference matrix corresponding to the penalty term. 
Here I x is the identity matrix restricted to x G Q s while vi and v 2 , are 
vectors with components 

(vi)fc = Xs{x k )Bi{x k - tt) (61) 
(v 2 ) fc = -X s (s fe ) J B 1 (27T - x k ) (62) 

In addition, di and d 2 are column vectors which approximate the derivatives 
of a vector u as 

u x (tt) w dfu (63) 
u x (2ir) « d^u (64) 

For instance, a centered difference approximation to the derivative u x (2it) 
would have (d 2 )jv-i — — (2/i) -1 , (d 2 )i = (2/i) _1 and (d 2 )fc = for k = and 
1 < k < N — 1. Lastly, since the support of -Bi(x) is restricted to x < 1, the 
function i?i(a; — 7r) = for x > 7T + 1. Hence, the numerical derivative of 
B\(x — 7r) at x = 27r is zero (or similarly with Bi(2n — x) at x = 0) 

d^vi = (65) 
dfv 2 = 0. (66) 



Combining the orthogonality conditions (65)-(66) with the fact that I x Vi ) 



.2 



vi >2 , implies that u = Vi )2 are eigenvectors with corresponding eigenvalues 

Ai = -r]~ l {l - df vi) (67) 
A 2 = -^(l-d^vz). (68) 

All other eigenvalues of B then lie in the space perpendicular to vi and v 2 



resulting in either A = or A = —r)~ l . The eigenvalues (67)-(68 ) are therefore 
directly a result of the modified penalty term and depend specifically on the 
component values of d 12 . As a result, the products dj 2 v 12 e (0, 1] depending 
on how one builds the numerical derivative vector di 2 . 
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As an example, taking a centered difference approximation to the deriva- 
tive u x (2tc) yields 

b^v 2 = ^((v 2 )i - (v 2 )*_i) (69) 

« 0.5. (70) 

The second line follows since (v2)jv-i — while (v 2 )i ~ h because the 
function B[(0) = I. 

In general, the product df 2 v 12 will be a weighted average of the deriva- 
tives of Bi(x) on the left and right of the interface. As a result, all eigenvalues 
A of B satisfy — 77 < A < 0. Therefore, modifying the penalty term does 
not change the time step restriction At < 2rj for a simple Euler scheme 

u n+l = u n + Af Bu n 

5.2. Numerical eigenvalues 



In the follow section we numerically compute the eigenvalues of (56)-(57) 
using a finite difference scheme for the spatial derivatives^ The scheme then 
has the form 



u n+1 = [l + At (L + B)]u n (71) 
where L is the standard 3-point stencil discrete Laplacian. As a result, the 



eigenvalues of the linear system ( 71 ) approach the real values associated with 
the Laplacian A when rj — > oo, and the values associated with the penalty 
term when 77 — >■ 0. 

To compute the eigenvalues numerically, we fix a grid with N points 
(256 < N < 4096) and examine the range 10~ 9 < 77 < 1. 



Proposition 2. (Practical stability) In practice, the numerical scheme (71) 
is stable provided one takes the time step restriction 

At < min{0.5/i 2 , 1.2^}. (72) 



Remark 6. The exact constant 1.2 in (|72j) depends on numerical details 
such as how one interpolates derivatives to the interface or the nature of the 
functions -Bo(^), Bi(x) 6 

Here figured shows that the numerical eigenvalues for iV = 2048 and rj = 10~ 7 



are stable with a time step restriction ( 72 ) 



3 Although not shown, a similar result of At < min{7V 2 , l.lrj} holds for a Fourier 
scheme. 
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Figure 4: Scatter plot of the numerical eigenvalues for (71). Here N = 2048, 
1] = 10~ 7 and At is taken from (72). 



6. Numerical example: Heat Equation 

In the following section we provide numerical examples for the heat equa- 
tion in dimension D = 1, 2. Specifically, we combine the analytic convergence 
and stability results from the previous sections to show that one may achieve 
high order numerical schemes. As a starting point, we demonstrate high or- 
der convergence in D = 1 dimension. We then move to D — 2, and outline 
additional details that arise from the numerical construction of the extension 

6.1. ID Heat Equation 

To test the convergence rates for the penalized heat equation, we use a 
manufactured solution approach. We note that the forced heat equation on 
x E [0,2vr], 

d t u = u xx + f (73) 



/ 

u(x, 0) 



sin(x+t) 
;,sin(x) 



[ cos(x + t) + sin(x + 1) — cos 2 (a; + t 



(74) 
(75) 



has an exact solution u e = e sm ( x+t \ To quantify the total error, we penalize 
d t u v = (u v ) xx + f - r/^Xs {u v - g) (76) 



equation (73) as 
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and compare the numerical solution of (76) to the exact one from (73)j 

To discretize in space, we use an equispaced grid with fourth order stencils 
for all derivatives. In addition, we treat all terms explicitly in time with a 
second order (improved) Euler scheme. When constructing the extension g, 
we first compute the derivatives of u at each grid point, i.e. u x (xk) or u xx {xk). 
We then interpolate the values of u x and u xx from the regular grid points to 
the points xr on the interface. 

Remark 7. The solution to the penalized heat equation u has a discontinu- 
ous second derivative u xx across the interface. As a result, interpolating u xx 
using regular grid points on both sides of the interface will produce a weighted 
average of right and left derivatives u xx (0~) and u xx (0 + ) in the construction 
of g. We note that in practice, such a procedure does not appear to alter the 
final numerical convergence rate. 6 

For our tests, we choose a solid region centered at 7r to be fl s = [tt— 0.7, n+ 
0.7]. To satisfy the stability restriction, we then take At = 0.2 h 2 , h = 2ir/N 
and slave rj = 5 At so that all parameter values are fixed by the number 
of grid points N. For each N = 2 k , with 6 < k < 12, we then numerically 



integrate (76) up to a final time of T = 1. We repeat the procedure using 
0, 1 and 2 derivatives of u in constructing the extension g and compare the 
numerical solution to the exact one (i.e. that of the unpenalized problem). 
Here figure [5] shows the convergence rates for matching different derivatives, 
while figures [6] and [7] show a typical solution and the corresponding error 
respectively. 

6.2. 2D Heat Equation 

In the following subsection, we outline the numerical details for aD = 2 
scheme. Here we work with an equispaced, regular grid with N x N points 
(64 < N < 512), and immerse the boundary T. The main difference when 
moving to higher dimensions is how one computes the extension g(x). To 
illustrate the construction, we refer to figure [8j To build g(x) we first compute 
all appropriate derivatives at each grid point (both on Q s and Q p ). For each 
grid point x£Sl s within a distance I of T, we compute £(x) as the orthogonal 
projection of x onto T and s(x) = ||£(x)— x|| 2 . Using a regular 9 point stencil, 



4 One can also restrict the forcing / = f(l — Xs) to the physical domain, and obtain 
similar results. 
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Figure 5: Plot of numerical errors \\u — w^numlloo f° r different values of N. 
The three curves correspond to building g using 0, 1, 2 derivatives of u and 
result with convergence rates of 1,2,3 respectively. 




Figure 6: Plot of the numerical solution u Vtnum (thin line) with the extension 
g (dashed line) for N = 2048. Here Q s = [ir — 0.7, tt + 0.7] is the solid domain 



we then perform a polynomial interpolation of all required derivatives from 
the grid points to Using the interpolated derivatives at £, one can then 
compute the normal derivatives of u required in equation (24) to construct 
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Figure 7: Plot of the total error in solving the penalized equations matching 
1 derivative of u in the extension. Here, iV = 2048, T = 1 is the integration 
time, while Q s = [n — 0.7, tt + 0.7] is the solid domain. 



g(x) at each grid-point inside £l s . Figure [9] illustrates a typical construction 
of g(x). 

Remark 8. For computational efficiency, one can precompute and store the 
values of £ as well as the appropriate coefficients required to extrapolate 
derivatives to the interface T. 6 

As an example in D = 2, we take the computational domain Q to be a 
periodic square with side length 2n. For the penalized domain Q s , we take a 
circle of radius r = 1/2 and center (x c ,y c ) = (tt,tt). The physical domain is 
then Q p = Q\Q S . To perform convergence tests, we again use a manufactured 
solution where u e = [e sm( - x ^ + cos(y)] cos(t). Here we perform a convergence 
test for the penalty parameter 77. To compute the convergence rate, we fix 
N = 512 and vary 5 x 10~ 5 < rj < 10 _1 , so that discrete numerical errors 
are smaller than the ^-dependent error obtained by introducing the penalty 
term. For different values of rj, we then integrate the penalty equation for a 



time T = 0.1 and compute the error. Figure 10 shows the L°° error between 



the penalized equation and the exact heat equation as a function of 77. 
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Figure 8: Regular grid with interpolation points. 



4 




Figure 9: A sample 2D plot of g matching 2 normal derivatives of it(x). The 
plot is taken at t = for the heat equation tests. 
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Figure 10: Converence plot of the L°° error at T = 0.1 for the heat equa- 
tion tests. The plot shows curves for matching derivatives (triangles), 1 
derivative (squares), and 2 derivatives (circles). The straight lines compare 
the expected convergence rates of 0(^ 5 ), 0{rj) and 0(r] 15 ) respectively. 

7. Numerical example: 2D incompressible Navier-Stokes 

The primary difficulty when transitioning from a penalized heat equation 
to the penalized incompressible Navier-Stokes equations is the addition of 
the velocity divergence constraint. Other differences, such as moving from 
a scalar to a vector equation, or adding a nonlinear convective term do not 
pose new additional challenges to the penalized equations. Intuitively, the 
difficulty with the divergence can be outlined as follows. For the penalized 
heat equation, the active penalty term forces the function u to closely track 
the extension function g. When moving to a set of vector equations, the 
velocity vector u,, will closely track the term g inside the penalty region f2 s . 
However, the component-wise construction of g will in general be such that 
Vg 7^ 0. Consequently, to remain consistent, one should not force V-u,, = 
inside fl s but rather allow V • to loosely track V ■ g. 

The simplest approach for handling the divergence constraint is to replace 
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V • u,, = with a Poisson equation [T3J [EH [2H 125] for the pressure. 



d t u v = F — Vp — -Xs(x)(u r) — 
F = — u ■ Vu,, + /iAu, + f 
Here the Poisson equation is 



Ap 




for x G f2 f 
for x G fh 



(77) 
(78) 

(79) 



Note that V • F is singular on T. As a result, we take the right hand side 



of equation (79) to be V • F inside the open set Q p , with an extension to 



on f2 s . Since the domain f2 is doubly periodic, the equation (79) does not 



require boundary conditions. Furthermore, we also remark that equations 



(77)-(79) recover the divergence constraint V-u^ = in f2 p . For any x G Vl p , 



we can take the divergence of (77) and substitute (79) to obtain 



X G Qr 



(80) 



Therefore, provided the initial data satisfies V ■ u = 0, we have V ■ u,^ = 
for all time. 

7.1. Discretization in time 



Here we outline a pseudo-spectral scheme for solving the equations (77) 



(79). We note here that given a velocity field u^, the equation (79) explicitly 
defines the pressure as a function of u^, p = p[u v ]. As a result, the right 



hand side of equations (77)-(78) can be discretized using any time stepping 



scheme. For instance, treating each term explicitly, one may use a simple 
Euler or Runge-Kutta scheme. 

For a second order scheme in h, we take a first order discretization in time 



with a time step restriction of the form outlined in (72). Since the domain is 



27r periodic, we can use the Fourier transform to invert the Poisson equation. 

Algorithm 1. (First order in time, second order in space) 
1. Compute an intermediate velocity u™ +1 

1 



r) r) 



At 



V 



X*(x)«- 



(81) 
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2. Compute the pressure 



Ap 



n+l 



1 

At 



(v-g; +1 )(i 



Xs 



(82) 



Using the Fourier transform, p n (k) = J-'fp' 1 ] with k = (k x ,k y ) and 
k = |k|, one may solve the Poisson equation efficiently as 



p" +i (k) 



f 



k 2 At 



^[(V-< +1 )(l- Xs )]. 



(83) 



Note that we have written (83) as the product of (1 — x) an d V • at 



the discrete level and refer to equation ( 79 ) as the definition of p. 
3. Update the velocity 



u 



n+l 



< +1 - (At^f-zkp 



n+ll 



(84) 



With the exception of inverting the Poisson equation, all derivatives in algo- 
rithm are computed using finite differences. 

Remark 9. In order to guarantee second order spatial accuracy in algorithm 
([I]), g should match at least 1 derivative of u„. 4k 

To test the order of accuracy of the active penalty method, we again use 
a manufactured solution of the form u e = (u e ,v e ) and p e where 



u e = cos(x) sin(y) cos(t) 
v e — — sin(a;) cos(y) cos(t) 
p e = sin(2a;) cos(y) cos(t). 



(85) 
(86) 
(87) 



Given initial data corresponding to the exact solution, we numerically 
evolve the velocity and pressure p using the pseudo-spectral method out- 
lined in algorithm [Tj Here we match 1 derivative of in the construction 
of g. Figure 11 shows second order convergence of the velocity field in the 
L°° norm, as expected. Meanwhile, the pressure and the divergence converge 



are one order less. As an example, figures 12a- 12b show the typical error for 



velocity and pressure while 13a and 13b show the horizontal velocity field 
along with the horizontal component of the extension g ■ x. 
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Figure 11: Navier-Stokes convergence plot. Second order convergence in the 
velocity field (squares), while the pressure (triangles) and velocity divergence 
(circles) converge at first order. 

8. Conclusion 

In this paper, we outline how to construct high order penalty methods. 
We do so by first introducing an active penalty term for the heat equa- 
tion. When we increase the number of matched derivatives, we show that 
the penalty term improves the analytic convergence rate in terms of the 
penalty parameter. Secondly, we examine the numerical stability of the ac- 
tive penalty term. We show that it does not introduce additional stiffness into 
the equations or additional length scales that would need to be resolved. The 
combination of the high order convergence in the penalty parameter along 
with the numerical stability then leads to higher order numerical schemes. 
Lastly, we extend the penalized term from the heat equation to the incom- 
pressible Navier-Stokes equations. In particular, we show how to handle the 
divergence constraint on the velocity field. 

Although we have outlined a high order approach, there are still remain- 
ing issues that limit the practical feasibility of the method. For instance, 
at no point do we improve the smoothness of the solution u. In fact the 
second derivatives of u remain discontinuous across the curve T, although 
matching more derivatives in the active penalty term reduce the size of the 
discontinuity. As a result, Fourier methods still have a slow decay in the 
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Figure 12: 
T = 1. 



(a) Velocity error (b) Pressure error 

Error fields with N = 128 for the velocity and pressure after 




(b) The extension g ■ x 

Figure 13: The numerical velocity field for horizontal component u v ^ um 
along with the extension function. Here N = 128 and T = 1. 



(a) Velocity u v jv- 
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Fourier modes thereby limiting the ability to spectrally compute derivatives. 
In addition, interpolation of high order derivatives in the construction of g 
should be one-sided (i.e. from Q p ) while in practice one would prefer to use 
points on both sides of T. As a result, ongoing research includes improving 
the global smoothness of u while retaining the high order convergence. 
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